Efficient Monte Carlo methods for simulating 
diffusion-reaction processes in complex systems 

Denis Grebenkov 1 ' * 

1 Laboratoire de Physique de la Matiere Condensee, 

CNRS - Ecole Poly technique, 91128 Palaiseau, France 

denis.grebenkov@polytechnique.edu 

(Dated: May 1, 2013) 

C^) ' We briefly review the principles, mathematical bases, numerical shortcuts and applications of fast 

random walk (FRW) algorithms. This Monte Carlo technique allows one to simulate individual 
trajectories of diffusing particles in order to study various probabilistic characteristics (harmonic 
measure, first passage/exit time distribution, reaction rates, search times and strategies, etc.) and to 
5_l , solve the related partial differential equations. The adaptive character and flexibility of FRWs make 

them particularly efficient for simulating diffusive processes in porous, multiscale, heterogeneous, 

■^T ' disordered or irregularly-shaped media. 
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I. INTRODUCTION 



Diffusion in complex systems often invokes numerical simulations. Except few simple shapes (such as a disk or a 
sphere) for which diffusion equations possess explicit solutions [1, 2], one needs to resort to numerical methods that 
can be roughly divided into two groups: 

O ■ 1. Finite differences, finite elements, boundary elements, etc. A domain and/or its boundary are dis- 

r/3 1 cretized with a regular or adaptive mesh. The original continuous problem is then replaced by a set of linear 

equations to be solved numerically. The solution is obtained at all mesh nodes at successive time moments. 

Since the accuracy and efficiency of these deterministic numerical schemes significantly rely on the discretization, 

mesh construction turns out to be the key issue and often a limiting factor, especially in three dimensions. 

Q-f 2. Monte Carlo simulations. A probabilistic interpretation of diffusion equations is employed [3-9] to represent 

the original continuous problem as the expectation of a functional of an appropriate stochastic process. Many 

,-H ■ random trajectories of this process are then generated and used to approximate the expectation and thus the 

^- 1 solution. Since there is no discretization, neither of the domain, nor of boundary conditions, Monte Carlo 

techniques are flexible and easy to implement, especially for studying diffusion in complex geometries [10-12]. 

o 

In basic Monte Carlo simulations, one fixes a time step 5 and approximates each trajectory r(i) by a sequence of t/6 
independent normally distributed random jumps along each coordinate, with mean zero and variance 2D6. Note that 
other jump distributions may be used, e.g., discrete displacements to neighboring sites of a lattice, ballistic displace- 
ments in random direction, etc. At each jump, one has to check whether the trajectory remains inside the confining 
C*~) ■ domain. If not, the jump has to be modified according to the imposed boundary condition. For instance, Neumann 
boundary condition is implemented by reflecting the trajectory into the domain, while the simulation is stopped for 
^- ■ Dirichlet boundary condition. Partial absorption/reflections, trapping, splitting and other local mechanisms can also 
be implemented. These fixed-time step simulations are easy to implement but are inefficient in hierarchical or multi- 
scale porous media. In fact, the time step <5 must be chosen so small to ensure that the average one-step displacement 
\2D6 is much smaller than the smallest geometrical feature of the medium. Since most particles are released inside 
large pores, a very large number of steps, t/6, may be required for simulating each trajectory. This is the major 
drawback of basic Monte Carlo schemes. 

In order to overcome this limitation, Miiller proposed the concept of variable time-step, geometry-adapted or fast 
random walks (FRWs) [13]. This technique was broadly employed by many authors, for instance, to simulate diffusion- 
limited growth phenomena [14-16] and to study diffusion-reaction processes and related first-passage problems in 
random packs of spheres [17-19] or near prefractal boundaries [20-23]. The idea is to replace Brownian motion by 
an equivalent "spherical process" that explores the medium as fast as possible. A particle starts to diffuse from an 
initial position Tq which may be prescribed (fixed) or chosen randomly. One draws the largest disk (or ball in three 



'Electronic address: denis.grebenkov@polytechnique.edu 



FIG. 1: A fast random walk in a medium with obstacles (dark disks). From an initial position ro, one determines the distance 
£o to the obstacles (or to the boundary of the medium) and draws a circle of radius £o centered at ro. A "jump" to a randomly 
(uniformly) chosen point on the circle is then executed. This single large displacement (shown by an arrow) replaces a detailed 
simulation of Brownian trajectory inside the disk of radius £q. From the new point n, one determines the distance again and 
executes the next jump, and so on (only three jumps are shown) [24]. 

dimensions) which is centered at ro and inscribed in the confining medium (Fig. 1). Its radius £q is the distance 
between ro and the boundary. After wandering inside the disk during a random time t\ , the particle exits the disk at 
random point ri. Since there was no "obstacles" inside the disk, all the exit points of the disk are equally accessible 
for isotropic Brownian motion so that the exit point ri has a uniform distribution on the circle of radius £q. From 
ri, the new largest disk of radius l\ is inscribed in the medium. After wandering inside the disk during a random 
time T2, the particle exits at random point r 2 , and so on. Following the Brownian trajectory of the particle, one can 
construct the sequence of inscribed disks (i.e., their centers r„ and radii £ n ) and the associated exit times r„. 

The fundamental idea behind FRWs is that the sequence {r„, £ ni r„} can be constructed directly, without simulating 
the underlying Brownian trajectory at all. At each step, one determines the distance £ n between the current position 
r„ and the boundary of the medium and chooses the next position r„+i randomly and uniformly on the circle of 
radius £ n . The time T n +x needed to exit from the disk (i.e., to jump from r„ to r„ +1 ) is a random variable which 
can be generated from the well-known probability distribution (Section II). A detailed time-consuming simulation 
of a Brownian trajectory with high spatial resolution is therefore replaced by generation of random jumps which 
are adapted to the local geometrical structure of the medium. In other words, the spatial resolution of the spherical 
process is constantly adapted to the distance to the boundary: closer the particle to the boundary, finer the simulation 
scale. Performing each jump at largest possible distance yields a tremendous gain in computational time. 

In this chapter, we review some practical aspects for an efficient implementation of FRWs and discuss several 
extensions and applications to first-passage phenomena. Section II introduces the formulas for generating first exit 
times and positions. In Section III, several strategies for estimating the distance to the boundary are discussed. 
Finally, Section IV presents extensions, applications and conclusions. 

II. EXIT TIME AND POSITION DISTRIBUTIONS 

We start by recalling the computation of the first exit time and position distributions for a general bounded domain 
£1 C M. d . Since the derived distributions will be used to generate jumps in a FRW algorithm, the underlying "jump" 
domain Q has to be simple (e.g., to be a disk or a sphere, as described in Section I). We emphasize that the simple 
shape of Q and the specific condition on its boundary dQ have nothing to do with the shape of the medium and the 
boundary mechanism that are going to be modeled with the FRW algorithm. 

We introduce the diffusion propagator G((ro, r) (also known as heat kernel, or Green function of diffusion equation) 
as the probability density for Brownian motion to move from a point ro to a vicinity of a point r in time t, without 
hitting the boundary dQ of the domain Q during this time. The diffusion propagator satisfies the diffusion equation 

^G t (r„,r) = i}AG ( (r ,r), (1) 



where A = d 2 jdx\ + ... + d 2 /dx 2 , is the Laplace operator, and D the diffusion coefficient. This equation is completed 
by the initial condition with Dirac ^-distribution, G t= o(r ,r) = 5(tq — r), stating that r is the starting point, and 
by Dirichlet boundary condition, Gt(ro,r) = at r <G dft, that "excludes" Brownian trajectories that hit or crossed 
the boundary dfl prior to time t. The diffusion propagator can be expressed in terms of the (Z^-normalized) Laplace 
operator eigenfunctions u m (r) and eigenvalues A TO [1, 2, 26] 

oo 

G t (r , r) = Y, e- DAmt w ro (r X(r), (2) 
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where asterisk denotes the complex conjugate, and the eigenvalue equation Au m (r) + X m u m (r) — is completed by 
Dirichlet boundary condition u m (r) — at r G d£l. 

The diffusion propagator allows one to compute the first exit time and position distributions. For Brownian motion 
started from a point ro, we denote r = inf{£ > : r(i) ^ 51} the first exit time r from the domain $7. Since r is the 
first exit time, Brownian motion should not hit the boundary dfl at earlier times t < r. The probability of this event 
(so-called survival probability) is obtained by integrating the probability density Gt (ro , r) over all the arrival points 
r: 

F ro {r > t} = S ro (t) = jdrG t (r ,r). (3) 

h 

The probability density is simply 



p ro (t) = -^- t S ro (t) = -JdrDAG t (ro,r) = Jds (-D^-G t (r ,r)) . (4) 

n an 



where we first used the diffusion equation and then Green formula to integrate by parts. Here d/dn is the normal 
derivative directed outwards the domain $7. The probability density p To (t) characterizes the time of the exit, regardless 
its position. Alternatively, one can consider the exit position s € dfl, regardless its time, whose distribution is known as 
the harmonic measure [25] . The harmonic measure density, p ro (s) , can be expressed through the diffusion propagator 
as 



oo 

p ro ( s )=Jdt(^-D^-G t 



(ro,r) . (5) 



Looking at Eqs. (4, 5), one can see that the diffusive flux — D^G t (r ,r) plays the role of the joint probability 
density for the exit time and position, while p ro (t) and p To (s) are marginal densities for time and position obtained 
after integration over the other variable. When both the exit time and position are needed, one can use the joint 
distribution. In practice, one can first generate the exit time and then generate the exit position conditioned to be at 
the prescribed exit time t: 

>"- le) = 7ZW>{- D l GM )„ (6) 

where the prefactor l/p ro (t) reflects the conditional character of this distribution. Finally, one may also need the prob- 
ability density of the arrival positions r inside the domain conditioned to survive up to time t, namely, Gt(i"o, r)/S ro (t). 
The above probability densities can be analytically computed for simple domains for which the Laplace operator 
eigenfunctions are known explicitly. In the following three subsections, we summarize the formulas for interval, disk 
and sphere. For convenience, we will use the rescaled (dimensionless) time and space variables, t — > Dt/L 2 and 
x — > xL, where L is the length of the interval, or the radius of the disk/sphere. 

A. Interval 

For the unit interval (0,1), the eigenfunctions are eigenvalues of the Laplace operator with Dirichlet boundary 
condition are 

u m (x) — v2 sin(7rma;), A m = Tr 2 m 2 (m — 1, 2, 3, . . .), (7) 



so that the one-dimensional propagator is simply 

oo 

E2 2 , 
sin(7rmx ) sin(7rmx) e _7r m . (8) 

m— 1 

It is also useful to write an alternative representation by image method: 

1 oo 

G t (x ,x) = 7 = J2 [e- {x - xo+2k f^ - e -(-+-o+2fe) 2 /(4*)] _ (9) 

k— — oo 

The first expansion (8) rapidly converges for large t, while Eq. (9) rapidly converges for small t. 

Using Eq. (3), one gets two alternative representations for the cumulative distribution function for the exit time: 



S Xa (t) = 2 ^ sin(7rmx ) 



l-(-l)' 



— 7r m t 



', x ,u - (10) 

,. , k-Xo\ , I k + x a 

erfc — =^- — eric 



i-«*^j + D-')* 



'At ) \ VAt 

where erfc(x) is the complementary error function. The density is 

oo 

p X0 (t) ^2ttJ2 sin(7rmx )(l - (-l) m )m e~^ mH 



k_£:c-i>w»).p(-<2±£). 



V4^ k= _ x 

As the boundary of an interval consists of two endpoints, the harmonic measure density is reduced to two functions, 
namely, the probability p Xa (s = 0) to reach the left endpoint before the right one, and that of the complementary 
event: 

v—v sin(7rmxo) , ,.. . 

fa s = 0=2V J y- = l-x , p xo (s = l)=x . (12 

^— ' 7rm 

m—l 

This is the classical result for the gambler's ruin [6]. Other probability densities can also be explicitly written. 

B. Disk 

For the unit disk with Dirichlct boundary condition, the Laplace operator eigenvalues and eigenfunctions are [1, 2] 

Kk = a? nk , u nk {r,6) = -£= — — J n (a nk r) cos nd, (13) 

where e„ = \/2 for n > and eo = 1, J' n {z) is the derivative of the Bessel function J n (z) of the first kind, and 
{<^nfe}fc=o,i,2,... is the set of all positive zeros of the function J n (z) (with n = 0, 1, 2, . . .). For convenience, the double 
index nk is used instead of the single index m to enumerate the eigenfunctions and eigenvalues. The asymptotic 
behavior of zeros a n k is well known while their numerical computation is straightforward by bisection method or 
Newton's method. After this preliminary step, the spectral decomposition (2) yields the diffusion propagator G t (r , r) 
and the consequent probability densities: p ro (t), p ro (s) and p ro _t(s). Since the starting point ro is located in the center 
of the disk, only the terms with n — (invariant under rotation) do contribute. In particular, one gets 

OO _„2 f 

f^ aokJi(ttok) 

while the harmonic measure density is uniform: p c (s) = j-, as all the boundary points are equivalent from the center 
(here and throughout the text, the subscript V refers to the starting point at the center of the jump domain). 



C. Sphere 

The Laplace operator eigenvalues and eigenfunctions for the unit sphere with Dirichlet boundary condition arc [1, 2] 
x nk = al k , u nk i{r,e,<p) = — ^= — T j n {a nk r)P l n {cose)e ilv , (15) 

V27T -J n ( a nk) 



where j' n (z) is the derivative of the spherical Bessel function j n (z) of the first kind, Pn(z) the associated Legendre 
polynomial, and {a n k}k=o,i,2,... the set of all positive zeros of the function j n (z) (with n = 0, 1, 2, . . .). The last index 
I ranges from — n to n. The spectral decomposition (2) allows one to compute all the necessary distributions. When 
the starting point ro is located at the center, the formulas are simplified, e.g., 

oo 

5 c (t)-2^(-l) fe+1 e- T2fc2 * (16) 

fc=i 

(as ao/c = ft(k + 1), k = 0, 1,2,.. .), while the harmonic measure density is uniform: p c (s) = j-. Using the techniques 
of Laplace transforms and series summation, one can get an alternative representation for the survival probability as 

[27] 

S c (t) = l-^f]e- (2k -V 2 /^\ (17) 

v?rt k=1 

which rapidly converges for small t. 

D. Behavior of the exit time distribution 

By definition, the survival probability monotonously decreases from 1 at t = to as t goes to infinity, as illustrated 
on Fig. 2a. In the short-time limit, the exit probability 1 — S c (t) is extremely small, 

' 2e -V(i6t)^2(i _ st + 192i 2 + 0(t 3 )) (d = 1), 
I -S c (t)~{2e- 1 /^(l-t + 4t 2 + 0(t 3 )) (d = 2), (18) 

r 1/(tt) A ( d = 3 )> 

because very few particles can travel the distance from the origin to the boundary during a short time [86]. In turn, 
for large t, the survival probability decays exponentially, 

f|e-" 2 * (d=l), 

ScW* j^feye-^* (d = 2), (19) 

l2e- 2 * (d = 3), 

since it is unlikely for diffusing particles to avoid the encounter with the boundary during a long time. 

Figure 2b shows the probability density, p c (t) — Jji °f the exit time. This figure and the above asymptotic 

behaviors clearly indicate that the exit time is localized around its mean value which is equal to l/(2d) (and 1/8 
in ID). In particular, the probability that r does not belong to an interval (i m in,imax), can be made negligible by 
choosing i m j n and < max appropriately. For instance, if i m ; n = 0.001 and t ma x = 10, one has 

¥{t <£ (0.001,10)} = 1-S C (0.001) + S C (10) < 10~ 12d (d = 2,3). (20) 

As a consequence, the exit times beyond this interval can be ignored. 

E. Generation of the exit times 

The explicit form of Eqs. (14, 16) allows one to generate exit times by inversion method. Inverting numerically 
the function S c (t) (i.e., finding a function T c {x) such that S c (T c (x)) = x for any x between and 1), one obtains a 
mapping from random variables \ n with a uniform distribution on the unit interval, to the exit times 

r„ = §T C ( X „) (21) 
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FIG. 2: (a): The survival probability S c (t) and its short-time and long-time asymptotic behaviors (18, 19) for the unit disk. 
The dashed horizontal line at S c (t) = x illustrates the construction of the inverse function T c (x). (b): The probability density 
—dS c (t)/dt of the exit time and its short-time and long-time asymptotic behaviors derived from Eqs. (18, 19). The vertical 
dotted line shows the mean exit time 1/4 [24]. 

for the interval of length £ n , or for the disk or sphere of radius £ n , with a given diffusion coefficient D. The uniform 
random variables \ n are generated by a standard routine for pseudo-random numbers. 

In summary, two preliminary numerical procedures arc required for generating the exit times: 

1. In 2D, one needs to find a finite number of positive zeros {cuofc} of Bessel function Jo(z). The inequalities 
7rfc < aofc < 7r(fc + 1) allow one to search for a single zero on each interval (irk, irk + tt) by bisection or Newton's 
method. Since smaller times require larger truncation sizes, the asymptotic formula (18) can be used instead of 
the truncated series in Eqs. (14, 16) to improve the accuracy at short times. In one and three dimensions, this 
step is skipped because the eigenvalues are known explicitly. 

2. One constructs the function T c (x) as a numerical solution T c (x) = t of the equation S c (t) = x for a mesh of 
points x. The monotonous decrease of S c (t) ensures, for any x, the unique solution which can be computed by 
bisection method or Newton's method. 

Both procedures rely on classical numerical methods. Moreover, these procedures have to be performed once and 
forever while the stored values of the function T c (x) can be loaded before starting Monte Carlo simulations. As a 
consequence, the generation of the exit times during simulations is reduced, through Eq. (21), to a routine generation 
of uniformly distributed pseudo-random numbers. 



F. Boundary condition 



By construction, a generated sequence of disks or spheres may become arbitrarily close to the boundary of the 
medium but it never hits the boundary. It is therefore necessary to set a threshold distance e below which the 
trajectory is considered to hit the boundary (Fig. 1). This threshold is typically set to be much smaller than any 
relevant length scale. Since the average number of jumps needed to reach the boundary scales as ln(l/e), setting e to 
very small values would not significantly slow down simulations. 

After the hit, the next step depends on the boundary condition or mechanism to be modeled. The two simplest 
conditions, Dirichlet and Neumann, represent purely absorbing and purely reflecting boundaries, respectively. In 
the former case, once a particle arrives onto the boundary, it is immediately "killed" and the simulation of the 
trajectory is stopped. This absorbing condition may model various physical, chemical or biological mechanisms, e.g., 
relaxation of local magnetization on paramagnetic impurities dispersed on the interface in NMR experiments; one-way 
permeation through cellular membranes; chemical transformation on the catalytic surface, etc. [28-32]. Whatever 
the microscopic mechanism is, the interaction with the boundary changes the state of the particle and thus removes 
it from the transport process. In the opposite case of Neumann boundary condition, the particle is reflected back 
into the medium and continues its diffusive motion, without changing its state. One can also simulate the so-called 
Robin boundary condition for partially absorbing/reflecting boundary when absorption and reflection events are 
chosen randomly [22, 33, 34]. Note that the boundary may have variable properties in space and time, e.g., some 
parts of the boundary may be reflecting while the other absorbing (this is the typical situation for search problems 
when a target is located on the boundary [35, 36]); moreover, these properties may evolve with time, e.g., during 
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FIG. 3: When the particle has approached the reflecting boundary of the medium closer than a prescribed threshold e, it is 
"released" on a circle of radius R centered at the encounter boundary point r' n . If the released point A does not belong to the 
interior of the medium, one uses its mirror reflected point r n +i inside the medium. The radius R should be chosen as large as 
possible, but small enough in comparison to the characteristic length of the boundary (so that the boundary is almost flat at 
scale R) [24]. 

passivation/deactivation of catalysts [37] . One can also consider multi-compartment media with permeable boundary 
(e.g., living cells and the extracellular space). Once the trajectory hits this boundary, it may be reflected on both sides 
(with equal or different probabilities). Apart from these standard boundary mechanisms, one can model trapping of 
a particle for a fixed or randomly distributed waiting time [38], multistage surface kinetics [39], and surface transfer 
mechanisms [40, 41]. Finally, the FRW algorithm for bulk diffusion can be combined with other numerical techniques 
for modeling transport on the boundary. Such a flexibility of FRW algorithms makes them attractive for simulating 
various diffusive phenomena. 

The implementation of Dirichlct boundary condition is trivial: simulation of the trajectory is simply stopped after 
hitting the boundary. One can record the hitting position (for computing the harmonic measure [20, 21]), the hitting 
time (for computing the exit time distribution [23] ) , the traveled distance (for computing the time-dependent diffusion 
coefficient [42-46]), etc. 

Neumann boundary condition is much more delicate, as the trajectory has to be reflected back into the medium. 
In the simplest implementation, one just moves the current position to the interior point of the medium at a fixed 
(small) distance from the boundary. This reflection distance has to be larger than the hitting threshold e but much 
smaller than other relevant length scales. The choice of the reflection distance is a compromise between the accuracy 
and rapidity of simulations: a large distance would bias the results, while a small distance would lead to multiple 
small jumps after reflection. Another, much more efficient way to reflect back the trajectory is to make a relatively 
large jump from the boundary of the domain to a half-circle or half-sphere (Fig. 3). Indeed, if the boundary is flat 
(or may be approximated as flat up to a certain length scale), diffusion in the related half-disk (half-ball) with the 
reflecting base is equivalent to diffusion in the whole disk (ball). For instance, the distribution of the exit times is the 
same, while the distribution of the exit position is still uniform. An implementation of this reflection jump allows one 
to significantly speed up simulations. 

G. Incomplete jump 

If the reflecting condition is set over the whole boundary, the trajectory would never stop. One adds therefore a 
condition to stop the simulation when the time counter t = t\ + ti + . . . + r„ exceeds a prescribed time T. More 
generally, whatever the boundary condition is, one may need to stop the trajectory after or at time T. If the trajectory 
has to be stopped precisely at T (e.g., to compute the traveled distance at time T), the last jump has to be simulated 
differently. In fact, we are interested in the intermediate position r(T) of the trajectory which is started at the center 
of the jump domain at time to = t — T n < T and conditioned to exit from this domain at time t > T. The trajectory 
can be split into two independent parts: from to to T, and from T to t. The conditional probability density for the 
intermediate position r is 



G T -t (r ,r)p r (t-T) 

Pr (t-t ) 



(22) 
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FIG. 4: A FRW simulation of reflected Brownian motion inside a cross-shaped medium by using square jump domains (shown 
in red) whose centers (full circles) belong to the trajectory (shown by solid blue line and generated by basic Monte Carlo 
simulations). On the left, we show a zoom of the first square jump domain for which the distributions of the exit time and 
position can be explicitly found by solving the appropriate PDE equation with Dirichlet boundary condition (denoted by 'D'). 
Each jump starts from the center. On the right, we show the third jump domain for reflection from the boundary. The right 
vertical edge contains the starting point and coincides with the reflecting boundary of the medium. The mirror reflection across 
this vertical edge reduces the analysis to that on the square with Dirichlet boundary condition. 

where the first factor is the probability density for moving from ro to r during the time T — to (first part), the second 
factor is the probability density for exiting the domain at time t — T after starting at r (second part), while the 
denominator is the probability density for exiting at time t — to after starting at ro. The intermediate position r can 
be generated from this density. 

H. Extension: rectangular domains 

The boundary of many model or image-reconstructed structures is formed by horizontal and vertical segments 
(in 2D) or by rectangular plates oriented along three coordinate axes (in 3D). For such structures, it may be more 
convenient to split the trajectory into "elementary blocks" inside rectangular jump domains rather than spherical ones, 
as illustrated on Fig. 4 [47, 48]. The use of rectangles/parallelepipeds as jump domains allows one to significantly 
reduce the number of small jumps near the boundary of the medium. On one hand, the generation of jumps becomes 
more complicated because the arrival points over the edges are not uniformly distributed any more. One needs 
therefore to generate not only the exit times but also the conditional exit positions according to Eq. (6). On the 
other hand, Brownian motion inside rectangles/parallelepipeds is formed by independent one-dimensional Brownian 
motions along each coordinate, and the problem is essentially reduced to one-dimensional setting on the interval. In 
particular, the diffusion propagator G t (ro, r) in the jump domain Q = [0, L{\ x [0, Lq\ x ... x [0, Lj] C K d is factored 
out as 



1 
Gf°(r ,r) = J7 — Gnt/L^roj/L^Tj/Lj), 



(23) 



3=1 



where Gt(xo,x) is the one-dimensional propagator in the unit interval [0, 1] given by Eqs. (8, 9), with dimensionless 
time Dt/ Li 2 - and dimensionless space coordinates xo and x, representing the starting and arrival points. 
The cumulative distribution of the exit time is the survival probability in £7: 



Sg{t) = jdv Gf (r ,r) = f[S roj/Lj (Dt/L 
o. i= l 



(24) 



where S Xo (t) is given by Eq. (10). As usual, the probability density is p ro v) = — Ji'S'ro (£)• I* 1 tne case 0I the 
(hyper)cube (with Lj = L) with the starting point in the center (i.e., ro = [L, L, ..., L]/2), one gets 

S c (t) = [Si(Dt/L*)] d , Pc {t) ^^[S^Dt/L^-'p^Dt/L 2 ). (25) 



Once the exit time is generated, one needs to generate the exit position at this time. For the sake of clarity, we only 
consider the starting point from the center of the domain. Starting from Eq. (6) and skipping technical derivations, 
we obtain the conditional probability density of the exit point s to lie on the face which is orthogonal to the coordinate 
axis j: 

qj (t) * G Dt/Ll (l/2,s k /L k ) 

c,t() ~~ii ^wm/Li) ' (26) 

where the last factor is the product of conditional probability densities for (independent) displacements in the orthog- 
onal directions k ^ j, and 

v^d r-2 P^ Ut l L k> 
l^k=\ L k Si(Dt/L 2 k ) 

is the conditional probability to exit through one of two faces at direction j. When all lengths are equal to each other, 
L\ = ... = Ld = L, one gets qj(t) = 1/d (all faces are equivalent from the center). 

The displacements along each coordinate can be generated by numerically inverting the cumulative distribution 

X 

Ft{x) 



(28) 
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A technical difficulty is that this function depends on two variables x and t, i.e., one needs to invert this function 
versus x for various values of t as a parameter. Other methods such as rejection sampling or adaptive rejection 
sampling can also be applied to generate these random variables. 

As discussed in Section II F, the boundary condition has to be implemented. When the particle is close to the 
reflecting boundary, one needs to perform a reflection jump, as discussed in Section II F. In this case, one or several 
edges (faces) of the jump domain may be reflecting (Fig. 4). In order to generate the exit time and position, one can 
consider diffusion in the extended jump domain which is composed of the original jump domain and its mirror image 
with respect to the reflecting boundary. The previous formulas can be directly applied. 

III. DISTANCE COMPUTATION 

The adaptive splitting of the trajectory into independent "elementary blocks" relies on estimating the distance 
from any interior point to the boundary. We emphasize that the jump distance may in practice be smaller than the 
distance to the boundary, but it must not exceed this distance. In other words, the problem of computing the exact 
distance can be relaxed to a simpler problem of finding a lower bound of the distance. Depending on the studied 
geometry, this purely geometrical problem can be solved in different ways. In this section, wc briefly describe two 
strategies for efficient distance estimation. 

A. Self-similar (fractal) domains 

Fractals are often used as a paradigm of complex domains [49]. On one hand, fractals are very irregular shapes 
which exhibit geometrical details at various length scales, "exploding" the classical notions of length, surface area or 
volume. On the other hand, self-similar or self-affine hierarchical structures help to perform accurate theoretical and 
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FIG. 5: Three generations of the quadratic von Koch curve of fractal dimension In 5/ In 3 
replaces all linear segments by the rescaled generator (first generation). 



1.465. At each iteration, one 





FIG. 6: (a) Basic arrow-like cell which is divided into the rotated square and five small triangles. Once a Brownian particle 
arrived into the rotated square, the distance between its current position (full circle) and the boundary of the arrow-like cell (the 
generator) can be computed explicitly. Random jumps inside the rotated square can therefore be executed, until the particle 
either exits from the arrow-like cell, or enters into a small triangle, (b) In the latter case, the particle starts to "see" the 
geometrical details of the next generation. The rescaled arrow-like cell is then used to compute the distance to the boundary. 

numerical analysis on fractals. In particular, the self-similarity of von Koch curves and surfaces allows one for a rapid 
estimation of the distance from any interior point to these boundaries [20, 21]. 

To illustrate this idea, we consider the quadratic von Koch curve (Fig. 5). When a random walker is far from the 
boundary, it docs not "distinguish" its geometrical details. One can thus estimate the distance by considering the 
coarsest generation (Fig. 6). Getting closer and closer to the boundary, the random walker starts to "recognize" smaller 
and smaller geometrical details. But at the same time, when small details appear in view, the rest of the boundary 
becomes "invisible" . Consequently, one can explicitly determine the distance by examining only the local geometrical 
environment (see [21] for details). The advantage and eventual drawback of this geometry-adapted algorithm is the 
need for a specific implementation for each studied geometry. Relying on self-similarity of von Koch fractals, we were 
able to generate up to 10 10 random trajectories for highly irregular boundaries with up to 10-12 iterations (in 2D). 
Note that these boundaries present geometrical features at length scales from 1 to (1/3) 10 (i.e., over five orders of 
magnitude). 



B. Sphere packs 



Sphere packs are often used to model porous and granular media, colloidal solutions, gels and polymer networks. 
This is indeed a quite generic model as spheres (or disks in two dimensions) can be mono- or polydisperse, overlapping 
or not, impermeable or not, while their locations can be regular or random. From the numerical point of view, 
sphere packs are particularly attractive as the distance from any point to the surface of a sphere can be easily 
calculated by knowing only the location and radius of the sphere. When the number of spheres is large (say, above 
few thousands), the computation of the distance to the boundary of the pack as the minimum over all distances 
to individual spheres becomes too time-consuming. In that case, one needs additional "tools" in order to rapidly 
locate a limited number of spheres that are the closest to a given point (the current position of Brownian trajectory). 
Among such tools, one can mention "coarse-grained" distance maps [15], Whitney decomposition or Voronoi cells. 
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FIG. 7: Whitney decomposition into dyadic squares near a DLA aggregate (shown in blue), with a zoom on the right, 
finest resolution is set to the diameter of small disks forming the aggregate (courtesy by D. D. Nguyen). 



The 



We illustrate a Whitney decomposition which is a partition of the mediu into squares/cubes such that the size of each 
square/cube is proportional to the distance from that element to the boundary of the medium (Fig. 7). Once such 
a decomposition of the domain is constructed (prior to launching Monte Carlo simulations) , the distance estimation 
is a very rapid: one determines the square/cube to which a given point belongs, and takes the size of this element 
as an estimate for the distance. In practice, one uses dyadic division into squares/cubes, while the decomposition 
is stopped when a chosen finest resolution is reached. The adaptive character of Whitney decompositions or other 
"coarse-grained" distance maps allows one to rapidly compute distances for sphere packs with millions of spheres. The 
related FRW algorithm can be used for simulating growth phenomena such as diffusion-limited aggregation [14-16], 
diffusion-reaction processes, search problems, etc. 

IV. APPLICATIONS, EXTENSIONS AND PERSPECTIVES 

A. Diffusion-weighted magnetic resonance imaging 

Diffusion-weighted MRI is a widespread experimental technique in which the random trajectories of diffusing nuclei 
(e.g., water molecules) are encoded by applying an inhomogeneous magnetic field [50-52]. When the nuclei diffuse 
inside a heterogeneous medium, the statistics of random trajectories is affected by the presence of walls or obstacles. 
Although these microscopic restrictions are not visible at the spatial resolution of DMRI, these geometrical features 
are statistically aggregated into the macroscopic signal. Measuring the signal at different diffusion times and magnetic 
field gradients, one aims to infer the morphological structure of a sample and to characterize the dynamics of a system. 
The non-invasive character of DMRI made this technique the gold standard in material sciences, neurosciences and 
medicine, with numerous applications to mineral samples (e.g., sedimentary rocks in oil industry; soils in agriculture; 
concrete, cements and gypsum in building industry) and biological samples (e.g., brain, lungs, bones). 

Until nowadays, modeling DMRI experiments was mainly restricted to simple structures [53-56], while truly multi- 
scale three-dimensional porous media such as concretes or sedimentary rocks remained out of reach due to the lack of 
efficient numerical techniques. The efficiency and flexibility of FRW algorithms make them promising tools for mod- 
eling DMRI. An approximate implementation of the gradient encoding into FRW algorithms was recently proposed 
that opens new opportunities in the study of mineral and biological samples [24, 57]. 



B. Continuous-Time Random Walks 



Many physical, chemical and biological transport processes exhibit anomalous features, e.g., the subdiffusive scaling 
of the mean square displacement [58-60]. Examples are the motion of organelles, vesicles or tracers in living cells 
[61-66], animals searching for food [67], contaminant or pollution spreading [68, 69], etc. For instance, overcrowding 
in living cells and colloidal gels, or deep wells in the interaction energy landscape [58, 70], may result in a heavy- 
tailed distribution of waiting times between jumps. Such long stalling periods would yield subdiffusive dynamics that 
can be described by Continuous-Time Random Walk (CTRW) with diverging mean waiting time but finite-variance 
spatial displacements. In that frame, the diffusion equation (1) has to be replaced by fractional diffusion equation, in 
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which the history of trapping or caging is incorporated through the fractional Riemann-Liouville derivative in time. 
Since the spatial dynamics is still governed by the Laplace operator, only the temporal dependence in the spectral 
representation (2) has to be modified as 

G t (r ,r) = J2 E a (-DX m t)u m {r )u* m (r), (29) 

771—1 

where E a (z) is the Mittag-Leffler function [71]. Using this extension, one can recompute the exit time and position 
distributions in order to simulate restricted subdiffusion in multiscale or complex structures [72-74]. 

C. Intermittent processes 

Another important extension of FRW is related to diffusive processes with two successively alternating phases 
(e.g., active and passive transport of vesicles in living cells [75-78]). These so-called intermittent processes have been 
intensively studied during the last decade. In particular, the alternation of phases was shown to facilitate search 
processes and speed up biochemical kinetics (see the review [41] and references therein). For instance, in the case of 
surface-mediated diffusion in spherical domains, the mean time for finding a target was shown to be minimal when 
the phases of bulk and surface diffusion alternate at certain rate [79-82] . At the same time, the question of optimality 
for the mean search time remains open for porous media or irregularly-shaped domains such as catalysts or biological 
structures. The FRW algorithm can be adapted to simulate such intermittent diffusive processes. 

D. Conclusions and Perspectives 

In this chapter, we described the basic principles and several applications of fast random walk algorithms. The 
algorithm relics on adaptive splitting of the random trajectory inside a complex medium into independent "elementary 
blocks" in simple jump domains for which simulations can be performed much more efficiently. In other words, the 
algorithm eliminates the geometrical complexity of the system and reduces the original problem to the study of 
diffusion inside a disk, a sphere or a rectangle, for which many mathematical tools (PDE, spectral methods, etc.) 
are particularly efficient [1, 2]. This "trick" drastically speeds up Monte Carlo simulations since at each step the 
largest possible exploration is performed. In practice, the computation is reduced to estimating the distance from any 
interior point to the boundary that can be solved in different ways depending on the studied geometry. The adaptive 
character of FRW techniques makes them well suited for simulating diffusion-reaction processes in complex, multi- 
scale, disordered, heterogeneous or irregularly-shaped media. The FRW algorithms can be applied for calculating the 
harmonic measures and their extensions, the first passage and exit time distributions, search times, reaction rates 
[83, 84], residence times [85], etc. Many growth processes such as diffusion-limited aggregation and related models can 
be efficiently realized by using these techniques. We also illustrated one application of these algorithms for computing 
DMRI signals and discussed extensions for modeling restricted anomalous diffusion and intermittent processes. 
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